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O ' Abstract 



In this note we apply some theoretical predictions that arise in the mean field frame- 
work for a large class of infinite range models to structural glasses and we present a first 
comparison of these predictions with numerical results. 



1 Introduction 

In recent years some theoretical progresses in our understanding of glasses have been achieved 
by comparing the results obtained for soluble models of generalized spin glasses |jl|, § with 
structural glasses properties under the assumption that the phase space of the two systems is 
similar. These progresses are also due to the use of the replica method and of related concepts 
as replica symmetry breaking, coupled replicas, dynamic transitions . . . 

In the nutshell many of the ideas of this replica approach are already present in the original 
papers of Gibbs and Di Marzio |^. However, the use of the whole panoply of tools developed 
in the study of spin glass models allow us to put these ideas in a much sharper form and to test 
them in numerical (and possibly real) experiments. Moreover the comparison of structural 
glasses and generalized spin glasses, introduced in ref . , has proven extremely useful and the 
conjecture that the two models have similar energy landscapes has been a fruitful starting 
point. 

In this note we will not discuss the theoretical basis under which this scenario has been 
derived but we will concentrate the attention on the physical picture and on the consequences 
of these predictions on the statistical properties of relative small samples. 

The basic assumption is that, at sufficiently low temperatures but still in the liquid phase, 
the system is almost always trapped for a long time in one of the exponentially large number 
of local minima of the free energy. The number (A/") of these local minima is related to the 
configurational entropy or complexity E by 

7Vr^exp(iVS(T)), (1) 

N being the number of particles and the previous formula being asymptotic for large A^. The 
total entropy density S is the sum of two contributions: the entropy density of each minimum 
and the complexity. This description is valid for Tx < T < T^. The complexity (which 
starts from a non zero value at To) is supposed to vanish linearly at a lower temperature, 
(i.e. at the temperature Tk) where the height of the typical barriers becomes infinite. The 
correlation time diverges at T^, and one can argue in favour of a Vogel-Fulcher law (e.g. 
r oc exp(— y4(T — Tk)~^), where v ^ 1). 

This scenario is implemented in a large class of infinite range models [§], where a more 
detailed picture of the phase space of the system is obtained. Indeed these infinite range 
models are soluble because the appropriate mean field theory is exact. The study of the 
mean field theory for these models is not a trivial task: the phase space structure of the 
configurations with low energy is quite complex and this is shown up in a rather interesting 
behaviour of the system both at equilibrium and in its approach to equilibrium. The existence 
of these complex structures implies the need of using modern tools (like the replica formalism) 
to study the properties of these systems. 

The aim of this note is to spell out some of the theoretical predictions that are obtained in 
the mean field framework for these infinite range models, to apply them (with the appropriate 
modifications) to the case of structural glasses and to present a first comparison of these 
predictions with numerical simulations. 

2 A mini theoretical review 

The main hypothesis of the Gibbs-Di Marzio's approach is that at low temperatures the system 
quite always stays close to a minimum of the free energy and that the long time dynamics is 



dominated by the time needed to escape from one valley to another. 

As noticed in ref. this scenario is implemented in long range models in which detailed 
computations can be done and a more precise picture can be obtained 0, |^. 

Let us consider for definiteness a spin model where the local variables are spins {cr{i), 
i = 1 . . . N). They could be either Ising spins a{i) = ±1 or spherical spins, i.e. real variables 
which satisfy the constraint 

E£^ ^ 1. (2, 

The Hamiltonian H{a) has a form that we do not need to specify here. There are many 
different models with quite different Hamiltonians which share the properties that we are 
going to describe. In the following we will consider a Gnite large system and only at the end 
we will send the value of A^ to cxd. 

We suppose that the phase space can be broken into many valleys separated by high 
mountains. In other words, the free energy (or the energy at low temperature) has many 
minima and the free energy barrier which we have to cross for going from one state to another 
is quite large. In the infinite range approximation it will be proportional to N. We can consider 
all the local minima of the energy as functions of the configuration and we can associate a 
valley to each of them at zero temperature. At higher temperature we have to consider the 
minima of the free energy and when increasing the temperature the valleys will disappear. 

We could also define the valleys in a dynamical way as regions of configuration space in 
which the system remains trapped for a large time. This definition is similar to the previous one 
if we assume that the large time behaviour of the dynamics at low temperature is dominated 
by activated processes which correspond to the crossing of high free energy barriers. 

Each valley (which we label by a) may be characterized by a magnetization 

ma{i) = (o-(^))a, (3) 

where {■)a is the statistical expectation value restricted to the a- valley. One can define a 
free energy F[m], which is a function of all the local magnetizations, whose functional form 
depends on the system. The valleys are local minima of this free energy. 

2.1 Equilibrium properties 

One finds in a very large class of models the following behaviour for the equilibrium properties 
depending on the temperature 0- |T0[] : 



• For temperatures higher than T^ there is only one minimum of the free energy F[m] at 
m{i) = Vz. In this case the total free energy of the system J-' is given by JF = F[0]. 

• At temperatures lower than Tm there is an exponentially large number of minima. The 
contribution of these minima to the partition function can be estimated by 

Z^ = Jdfexp{-Nf3f + ArS(/,/3)) ^ exp^-NPf* + iVS(r,/5)), (4) 

where / is the free energy density {Nf = J-'), exp(A^S(/, /?)) is the number of minima of 
the free energy density / and /* is the value of / which maximizes the exponent, that 
is a function of /3. 

Below Tm we can distinguish three regions: 



• For T > To the contribution to the partition function of the non trivial minima can be 
neglected and the free energy is still given by F[0]. 

• For To > T > Tk the contribution to the partition function of the non trivial minima is 
dominant. The number of minima which dominate the partition function is exponentially 
large and the total entropy of the system is given by 

S = S^ + NE{f\(3), (5) 

where Sm is the contribution to the total entropy of one minimum and S(/*,/3) > 0. 

The magnetization averaged over all the minima is given by m{i) = (cr(i)) = J2a Wa{<y{i))ci 
where Wa oc exp(— /JF^j) is proportional to the contribution of the a-minimum to the 
partition function (SqWq, = 1). 

In this region the magnetization averaged over all the minima is zero and the total free 
energy which is dominated by the contribution of all the non trivial minima is still given 
by -F[0]. The free energy and all the other static equilibrium properties (quantities?) 
are fully regular at To- 

• For T < Tk, W*,P) = 0- For T slightly larger than Tk one finds S(/*,/?) oc (T - 
Tk), so that the complexity (and consequently the entropy) has a discontinuity in the 
temperature derivative at T = Tk, which is the transition point from a thermodynamic 
point of view. Here the entropy becomes equal to the contribution of a single minimum 
and this temperature may be identified with the Kauzmann temperature. 

The region T < Tk can be characterized by a peculiar behaviour. The partition function 
is dominated by those states which have the minimum free energy and there are some minima 
for which quantities Wa remain of order 1. 

In the previous description we have neglected the possibility of crystallization, i.e. the 
formation of a highly ordered state which leads to a first order transition. If we consider this 
possibility we have to distinguish two cases, systems with quenched disorder in the Hamiltonian 
and systems without disorder (often with a translation invariant Hamiltonian). 

• There are systems whose Hamiltonian contains quenched disorder. A typical example 
would be: 

H = Y,J{hk,l)a{i)a{k)a{l), (6) 

i,k,l 

where the variables J are random (e.g. Gaussian). 

Another example (in which the mean field approximation is non exact) would be a system 
in which the particles interact not only with themselves but also with an external, fixed, 
random potential. It is clear that if the external potential is strong enough the free 
energy of the crystal phase may become quite large, while that of the glassy phase may 
be much less affected. 

• Other systems may not contain quenched variables in the Hamiltonian. A typical exam- 
is: 
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H = ^Y: sin [-j^^ ia{^) - a{k)f with a(.) = ±1. (7) 



Another example would be a system in which the particles interact only with themselves 
with a given potential. 

The two categories of models seem rather different one from another; however it has been 
noticed that in the mean field approximation systems belonging to the two categories behave 
in a quite similar way, with the only difference that systems without quenched disorder may 
crystallize. 

While this kind of results can be proved in the long range models, their validity for short 
range models (like structural glasses) remains an open question. In this note we would like to 
present numerical evidence for the correctness of these ideas also for structural glasses. Before 
doing so we will examine in detail the predictions of the mean field approach. 

2.2 Equilibrium properties at low temperature 

For simplicity let us consider the predictions that would follow from the application of the 
replica approach to a system of A^ particles interacting with a given Hamiltonian in a box. 

We will consider for simplicity a system without any type of symmetry (the box is not 
symmetric and the particles interact with the walls, periodic boundary conditions are not 
used so that there is no translational invariance). We suppose that for a given value of A^ the 
Hamiltonian has many minima which we label by a. 

The partition function can thus written as 

Z(/3,Ar)^5:exp(-/5F«(AAr)) 
F„(/5,A^)=E„,^-T5(a,Ar), (8) 

where for simplicity we label the configurations in increasing free energy order. For small 
temperatures the entropic contribution to the free energy will be as usual negligible. It is 
quite evident that for a finite system the partition function is dominated by the lowest energy 
configurations in the limit of zero temperature, and we will suppose that this property persists 
in the infinite volume limit. 

For different values of A^ we will have quite different values of the free energies. The lowest 
free energy states will have free energy differences of order 1, so that by adding a single particle 
(i.e. going from A^ to A^+ 1) we can strongly change the values of the Fa{(3, N). If we label by 
Fq the lowest free energy we expect that the quantities Fa{P, N) — Fq{P, N) remain of order 
1 when A^ goes to infinity, but they do not go to any limit (they have no A^ — *> oo limit?) 
because they change with A^. In this case it is natural to introduce a probability distribution 
for the differences in free energies, which tell us the probability of finding a given value of 
Fa{f3, N) — Fq{(3, N) if we chose randomly (albeit large) the value of A^. In other words if the 
value of the free energy changes strongly with A^ it is convenient to introduce the probability 
distribution of their values. 

The precise construction is the following: for any large value of A^ we introduce a reference 
free energy Fr(/3, A^) such as that, in the whole region F — Fp>{i3, N) << N, the probability 
of finding a minimum of free energy F is given by 

PiF) = exp{P{x{P){F-Fn{P,N))), (9) 

where the quantity x{j3) parameterizes the probability distribution. 



The condition that the free energy is approximately given by Fr imphes that the integral 

/ P(F)exp(-/5(F-F^(AiV)) (10) 

is convergent and therefore x{f3) < 1. 

The probability of finding a configuration in the minimum labeled by a is 

Wc, oc exp(-/5F«), (11) 

where the normalization condition 

E^« = l (12) 

a 

is satisfied. 

In many models x becomes equal to 1 at the temperature Tk- At low temperatures x is 
proportional to the temperature (it would be exactly proportional to the temperature if we 
neglect the entropic contribution). 

An important property of the model is how different the various minima are. The simplest 
hypothesis is that they are as different as possible, i.e. the correlations among the particles in 
one minimum and in another minimum (among those of lowest energy) are zero. For example 
the probability of finding two particles in two different minima at a given distance r does not 
depend on the distance: 

-Pa,7(r) = ^^ {xa{i) - Xy{k) -r) = p^ (13) 

i,k 

This hypothesis is usually called one step replica symmetry breaking. More complicated 
distributions of the distances are discussed in the literature. 

If there are symmetries the situation becomes slightly more complex. For example in 
a translational invariant system if x{i) are the coordinates of a minimum, x{i) + 6 are the 
coordinates of another minimum. It is therefore useful to consider all the minima which are 
related by a symmetry transformation as a single minimum. 

2.3 The probability distribution of the distance 

The properties of the system may be sharpened by introducing an appropriately defined dis- 
tance d between configurations of particles and looking at the corresponding equilibrium prob- 
ability distribution: 

PN{d) = X! U]aU!^6{d — da-y)- (14) 

a, 7 

In the limit of large A^ we have that 

PN{d) = a^Sido — d) + h]\[5{di — d), (15) 

where a^v + ^at = 1- Of course for a finite system the delta function will be smoothed. 

The function a^ represents the probability of finding two different configurations in the 
same minimum and it is given by 



In the one step replica symmetry breaking hypothesis two configurations that are not in 
the same minimum are expected to be orthogonal, i.e. at the maximum possible distance di. 
This happens with the probability: 

bN =Y1 WaW^. (17) 

The functions oat and hjsf are naturally depending on the temperature. In the mean field 
framework one finds that the weight of the (5-function in di, h^., that is zero for T > T^-, 
increases continuously when lowering the temperature below the transition point. 

The probability distribution PN{d) should therefore be an appropriate observable for look- 
ing at the transition from the high-T region (where Pnid) is Gaussian-like for a finite system) 
to the glassy phase in which it is expected to show a non trivial behaviour, strongly depending 
on A^. 

2.4 The approach to equihbrium 

For simplicity we will point the attention on the relaxation of the energy density when the 
system is quenched abruptly (at the time t = 0) from a random initial configuration (i.e. 
infinite cooling rate) to the final temperature. 



In the mean field picture |T3[ one finds two different relaxation behaviours, depending on 
the final temperature value: 

• At high temperatures the energy reaches exponentially its equilibrium value, i.e. e(t) = 
Ceq + c exp(— t/r). The relaxation time r increases when lowering the temperature and 
it diverges for T ^ Td- 

• Below the dynamical transition point To the energy behaves linearly as a function of t~" 



1I4| , i.e. e{t) = eD + ct'"", the exponent a being weakly depending on T. Here the system 
is evolving towards some metastable states (with infinite life time) and correspondingly 
the asymptotic energy value cd is higher than that of the equilibrium. 

The infinite life time of metastable states is just an artefact of the mean field approximation. 
In a real system we expect the approach to the equilibrium for T <Td consisting of two steps: 

• The convergence to some metastable states with a mean field like behaviour. 

• The slow decay of metastable states due to activated processes. In this second step the 
system reaches the true equilibrium state that will be still the replica-symmetric one for 
T>Tk. 

This means that by looking at e(t) it is in principle possible to find numerical evidence for 
the reminiscence in real glasses of the mean field dynamical transition. Here T^ is expected 
to mark the onset of the two steps relaxation. 

3 The model 

We study a binary mixture of soft spheres, half of the particles being of type A with radius 
a A and half of type B with radius as- The Hamiltonian is: 

n,.. £ (;i(M))". „8, 

i<k=i \ r^ ^^=1 / 



This model has been carefully studied in the past ||2D|-|2^. It is known that the choice 
<^b/<^a = 1-2 strongly inhibits crystallization. We also follow the convention of considering 
particles with average diameter 1 by setting: 

a\ + 2{aA + aBf + al = ]^. (19) 

Thermodynamic quantities only depend on F = p/T^/^^ where T = 1/(3 is the temperature 
and we have taken density p = 1 (F = 13^^^)- The A^ particles move in a ?>d cube of size 
L = N^^^^\ The glass transition is known [^] to happen around F^ = 1.45. 



In order to obtain numerical results on the equilibrium properties comparable with the 
mean field theoretical picture we attempt to measure the equilibrium probability distribution 
P{d) of the distance d between states. Following the usual strategy in spin glass simulations 
we introduce two evolving contemporaneously and independently replicas of the system. P{d) 
is then given by: 

P{d)^{5id~V)), (20) 

V being the appropriately defined distance between the configurations of the two replicas. 

Labeling by {vi}, {si} the positions of the A^ particles in the two replicas, a natural 
definition of T> is the Euclidean one, minimized over permutations tt, rotations and reflections 
R and, when using periodic boundary conditions, translations T: 

^' ^ ^ "^^^T I E (^^ - ^A^)y + E (^^ - ^A^)n ■ (21) 

iV TT.H.l ytyp^A typeB / 

An analogous definition of distance between configurations of particles has been considered in 



p!5|] and in |T6|, in studies of potential energy minima for Lennard- Jones systems. 

We minimize over permutations by an approximate procedure. For each particle i of one 
configuration we take iT{i) being the nearest one of the other. We expect this to be a reasonable 
approximation in the considered temperature range (from F=l to F=2) since the probability 
for two particles of the same system of being at a distance lower than 2(T ~ 1 is very close to 
zero at not too high temperatures (the radial density- density correlation function in a simple 
liquid shows the well known behaviour g{r) ^0 for r < 2cr). 

Minimization over rotations and reflections is easily performed since in the case of a 3d 
cube R results a discrete group that includes 48 symmetry operations, corresponding to the 
2^ reflections and to the 3! permutations of axes. On the other hand, minimization over the 
continuous group of translations is a hard task and we prefer to avoid it by not using periodic 
boundary conditions. 

To measure P{d) we have considered a slightly modified model in which particles are 
definitely confined in the cubic box of size L by a soft walls-repulsive potential term: 

TV 3 / ^ 1 \ 

We have chosen the values Ci = 7r/5 and C2 = 0.6 that give a behaviour of the energy as a 
function of F quite near to that of the periodic boundary conditions case. 



4 Numerical results on the dynamics 

4.1 Algorithms 

Both for finding the best numerical approach to the simulations of structural glasses at the 
equilibrium and for studying the dynamical properties of the model we have implemented 
different algorithms, considering stochastic and deterministic dynamics: 

• Monte Carlo {MC). We start from a random configuration and we quench abruptly the 
system by putting it at the final temperature (i.e. infinite cooling rate). During one 
step each particle is suggested to move of a random quantity and the maximum shift 
permitted is chosen in order to obtain an acceptance ratio near 0.5. 

• Molecular Dynamics (MD). When using MD to start with a completely random spa- 
tial configuration may cause difficulties. Here the initial spatial configuration is obtained 
from a random one by 60 Monte Carlo steps at the final temperature. At the begin- 
ning and again each 100 steps we extract momenta {pj} according to the Boltzmann 
distribution (we have taken the masses of particles ttia = itib = 1) and we impose the 
condition Y^iLiPi^ = (3iV — fMD)T, where /md is the number of frozen degrees of free- 
dom. Jmd = 3 when using periodic boundary conditions since the three components of 
the total linear momentum P are conserved (we have taken P = 0), otherwise Jmd = 0. 
We use the velocity- Ver let algorithm |l7] with St = 1/250 (a value that we have checked 
to be reasonable in our case, e.g. the total energy results perfectly constant). 

• Isothermal Molecular Dynamics (IMD). The system evolves according to the Gaussian 
isokinetic equations of motion |]T^ : 

V^ F rt 

h=P. P. = F,-Xp, ^= ^'=' ', ■ (23) 

2^i Pi 

At the beginning we extract momenta and we fix the constraint Y^iLi Pi^ = (3A^ — 
fiMD)T, where fiMD = fuD + 1- Also in this case we start from a spatial configuration 
obtained from a random one by 60 Monte Carlo steps at the final temperature. We 
implement the algorithm by using a variant of the leap-frog scheme [|l^ again with 
5t = 1/250. 

• Parallel Tempering {PT). Algorithms in which temperature is allowed to become a 
dynamical variable |]18| are very effective for thermalizing systems with a complex free 
energy landscape. In this recently introduced [|Ty| method a set of n different (3 values 
/?i < ■ ■ ■ < /3fc ■ ■ ■ /5„ is chosen a priori and n replicas of the system evolve contempora- 
neously. The extended Hamiltonian HpT = Sa=i Pi^) '^[Ca] is defined and exchanges 
of temperatures between replicas 

(/3(ai) = A„/3(a2) = /?fe) ^ {P{a,) = A„ /^(aa) = Pk,) 

are allowed with probability p = min[l,exp(— ATipr)]. The whole process is itself a 
Markov chain and, when equilibrium has been reached, each replica moves between 
different temperatures of the set remaining at the equilibrium. We start by extract- 
ing independently the n random initial configurations and we quench each replica at a 
different one of the chosen temperatures. In a PT step, sequentially for a = 1 . . .n: 

9 



— Replica a makes one MC step at its temperature (3{a) = (3k- 

— A random number j = ±1 with equal probability is extracted. 

— For 1 < k + j < n the exchange of temperatures between replica a and replica 6, 
where [3{h) = Pk+j, is suggested and possibly accepted. 

The set of (3 values should be chosen carefully. We find that for the numbers of particles 
considered (up to A^=36) PT works well down to r=2 with ?t,=13 different temperatures 
(r=l, 1.05 ...1.2, 1.3 ... 2). 

We have implemented and tested a combined technique of Parallel Tempering and 
Isothermal Molecular Dynamics (IMDPT). The n initial spatial configurations are 
obtained from random ones by 60 MC steps at r=l. After extracting momenta we 
impose the constraints Z]ili(Pi")^ = (3iV — fiMD)Ta a = 1 . . .n (we have taken the same 
set of temperatures that in PT) . Replicas evolve accordingly to IMD and each 10 IMD 
steps the exchanges of temperatures happen with probability p = min[l, exp(— AT-^p^)], 
where in HpT only potential energies appear. Here to change temperature Tdd — > T^ew 
means to transform momenta as p^ -^ {Tnew/ToidY^'^Pi, i = 1 . . . N. 

Finally we have performed MC simulations of large systems (2000 and 8000 particles) in 
the periodic boundary conditions case, by putting a cut-off on the soft spheres' potential: 






^) for \ri-rk\ < R 

Vik = { (24) 



( cr{i)+o- 



I = ^fc otherwise 



The algorithm is then implemented in such a way that for each particle the map of the 
ones which are at distance lower than i? + 25 is recorded and updated during the run 
{5 being the maximum shift permitted to a particle in one MC step). We have chosen 
R = 1.7 that means a practically negligible V^^ ~ O(10~^) {Vik ~ 1 for |rj — rk\ ~ 2a). 

4.2 On the energy relaxation behaviour 

We look at the (potential) energy density relaxation when the system is quenched abruptly 
from a random initial configuration to the final (low) temperature. We compare results ob- 
tained by different numerical methods for A^=34 particles and we extend the analysis to large 
systems (A^=2000 and iV=8000) in the MC case. Here periodic boundary conditions are used. 

In [Figs. 1-3] we present data obtained by MC, MD and IMD, respectively. Our results 
show no evident difference between stochastic and deterministic dynamics. Not only in the 
MC case, as already observed in previous simulations on the same model ||26|, but also when 
using Molecular Dynamics techniques the energy density relaxation is well compatible on a 
large time window with a linear behaviour in function of t~°. We get a non small exponent 
a ~ 0.8, weakly F-dependent in the considered range (the dependence seems slightly more 
pronounced in the IMD case). 

Moreover, the energy values at a fixed temperature that one can estimate by asymptotically 
extrapolating the linear behaviours obtained by different methods, are near to be the same. 
On the other hand we are going to show that these asymptotic values (apart from the case of 
r=1.4) are not the real equilibrium energy values of the system. 

10 
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Figure 1: MC data on e as a function oi y = t^""-^) at r=1.4(+), 1.6(x), 1.8(*) and 2.0(n) 
Here 150 different initial conditions are considered. 
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Figure 2: MD data on e as a function of y = t^""-^) at r=1.4(+), 1.6(x), 1.8(*) and 2.0(n). 
Here 100 different initial conditions are considered. When using Molecular Dynamics tech- 
niques the need of a non completely random initial configuration (i. e. the first 60 steps are 
MC steps) induces some short time effects and the linear behaviour takes place at a slightly 
larger time than in the MC case. 
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Figure 3: IMD data on e as a function oi y = t^""-^) at r=1.4(+), 1.6(x), 1.8(*) and 2.0(0). 
Here 100 different initial conditions are considered. The dependence of a on F seems slightly 
more pronounced in this case, our best estimates ranging from a ~ 0.8 at r=1.4 to a ~ 1.1 
at F=2.0. 
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We present in [Figs. 4-5] PT and IMDPT data. There is a reminiscence of the hnear 
behaviour but here the system moves between the high and the low temperatures of the 
set, this preventing it from being trapped. At the end of the (quite large) time window, 
equilibrium is near to be reached and the energy values are definitely smaller than the ones 
estimable by extrapolating simple Monte Carlo or Molecular Dynamics data. The difference 
is already detectable at r=1.6 and it becomes more pronounced at lower temperatures. To 
further outhne this result we plot in [Fig. 6] data obtained by various methods at r=2.0. 




Figure 4: PT data on e as a function oi y = t(-°-^) at r=1.4(+), 1.6(x), 1.8(*) and 2.0(0). 
Here 24 different initial conditions are considered. Note that both in this case and in the 
IMDPT one the obtained behaviour depends on the entire set of temperatures. 
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Figure 5: IMDPT data on e as a function of t(^°-^) at r=1.4(+), 1.6(x), 1.8(*) and 2.0(0). 
Here 24 different initial conditions are considered. Data are near to be compatible within the 
errors with the PT ones. 

We have found numerical evidence for the energy relaxation at low temperatures consisting 
of two processes that happen on remarkably well separated time scales. Both in the MC and 
in the Molecular Dynamics data the first step is clearly observable, corresponding to the 
convergence to some metastable states with a mean field like behaviour (note that a is weakly 
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Figure 6: Data on e as a function of t*^ °-^^ at r=2.0 obtained by MC{+), MD{x] 
PT(n) and IMDPT(o). 
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depending on the temperature on a large range). The fact that the extrapolated energy values 
are not the real equilibrium values gives evidence for the presence of a second step, the slow 
decay of metastable states dominated by activated processes. 

These results are inadequate for understanding whether the curvature of low temperature 
data on eiy) at very large times (small y = t~°) is related to the quite small N value we 
are considering or it actually represents the beginning of the second step (possibly with an 
intermediate plateau). In order to clarify this point we consider definitely larger systems. 

We plot in [Fig. 7] e{y) as obtained by MC simulations for iV=8000 at 1=1.4, 1.5, 1.8 and 
3.0. At r=1.8 we present data also for iV=2000 that result almost indistinguishable from the 
A^=8000 ones and remarkably similar to MC data for the well smaller A^=34 system. This 
means that we are actually looking both at the first relaxation step and at the beginning of 
the second one. 

The different behaviours observed when varying the temperature agree well with the pre- 
viously discussed theoretical picture: 

• At r=1.4 no two steps behaviour is observable and the energy reaches the equilibrium 
value in the considered time window. The system appears to be still above T^ but quite 
near to it (the relaxation time r seems to be very large, the expected exponential decay 
being not distinguishable from a linear behaviour in function of y). 

• For r > 1.5 the system is definitely below Tq. The energy relaxes linearly as a function 
of 1/ to a value cd higher than that of equilibrium, the slow decay of metastable states 
happening on a well larger time scale. 

• The beginning of the second step is clearly observable at r=1.5 and still quite evident at 
r=1.8 but it disappears at r=3.0. Here data behave linearly on the entire time window. 
This seems reasonable since the time scale that control the second step (i.e. the onset 
of activated processes) is expected to increase when lowering the temperature. 
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Figure 7: MC data on e as a function of t*^ '^•^-' respectively for A^=8000 (one initial condition) 
at r=1.4 (+), 1.5(x), 1.8(*) and 3.0(n) and for A^=2000 (two different initial conditions) at 
r=1.8 (o). 
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Figure 8: MC data on e as a function of t(-°-^) for A^=8000 at r=1.5 (+) and 1.8(x). The 
lines are the corresponding best fits to e(t) = at*^~°'^) + 6t''' + c obtained by the plotted points. 
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To fit the two steps behaviour is a hard task. In [Fig. 8] we show our best resuhs at r=1.5 
and 1.8 (A^=8000), obtained by fixing ot = —0.8 and considering e(t) = at^~°'^^ + bt"' + c It is 
interesting to note that in both cases we get a positive 7 value of the same order, 7 ~ +0.1. 
It seems reasonable that the second step corresponds to a streched exponential decay, i.e. 
e(t) = ai(t/ri)— + a2exp(-C(t/r2)+l^l) +eeg ^ at^-o-^) - \b\t+^^^ + Cpiateau for t < ra. However 
we stress that this is a purely indicative result since quite different 7 values (possibly of the 
opposite sign) are obtainable by slightly varying a or the time window. 

To conclude this section, we note that IMDPT works well but seems to be not more 
effective than usual PT. We have done some dynamical simulations with different methods 
also in the case of Hamiltonian (|22|). We do not show the results that are qualitatively similar 
to those obtained by using periodic boundary conditions, the only difference being that here 
the onset of the second step becomes evident at shorter times. 

5 Numerical results on the statics 

5.1 Simulations at the equilibrium 

In order to get the equilibrium behaviour of P{d) we use Parallel Tempering, simulating 
contemporaneously two independent sets of replicas. Thermalization is a hard task and we 
limit our analysis to quite small numbers of particles ranging from A^=28 to A^=36. For each 
A^ value we perform an extensive simulation of 2^^ -=- 2^^ PT steps (i. e. up to more than 16 
millions of MC steps for each one of the 26 replicas). We measure all the quantities we are 
interested in during the last 3/4 of the run. 
Thermalization is checked in different ways: 

• We divide the last part of the run, in which statistics is collected, into 16 equal intervals 
and we look for possible shifts of the corresponding mean values (particularly we do not 
find evident changes in the behaviour of P{d)). 

• We check that each replica moves more than once from an extrema of the temperature 
range to the other and back in the last part of the run. 

• We valuate the specific heath c both using c = d < e > /dT and using T^ c =< e^ > 
— < e >^, checking compatibility of results. 

Errors are estimated from the mean values obtained in each one of the 16 intervals. 

Despite our efforts it is not possible to exclude the presence of lower minima in the free 
energy landscape that are not accessible to the system on the considered time scale. To clarify 
this point it is interesting to look at the particular value A^=32=4^/2. Here particles are 
allowed to fill all the sites of a fee structure. Once this crystalline equilibrium state has 
been reached (in about 2^*^ -h 2^^ PT steps) the system rests definitely trapped in it, the 
crystalline configuration corresponding to a "golf-hole" shaped minimum in the free energy 
landscape. The energy density results a discontinuous function of F as expected for a first 
order (liquid-crystal) transition. 

It seems unlikely to us but particles could be able to arrange themselves in some sort of 
crystalline configuration that we have not detected, also for A^ 7^ 32. On the other hand we 
would like to look at the glass equilibrium states, apart from the possible presence of lower 
crystalline minima. We find therefore reasonable to assume that on the time scale studied the 
system extensively explores the entire phase-space region we are interested in. 
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5.2 On the behaviour of P{d) 

Data on the energy density e(r) for different N values are plotted in [Fig. 9]. We note that 
the equilibrium energy density is a not regular function of the number of particles since it is 
definitely smaller for iV=34 than for A^=36. The non trivial iV-dependence of the model at 
low temperatures is outlined when looking at the specific heath c(r) [Fig. 10]. 
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Figure 9: Data on the energy density e as a function of F for A^=28(+), 30(x), 34(*) and 
36(n). Lines are only to join neighbouring points. 

It should be stressed that for each given N value the energy appears a continuous function 
of F and correspondingly the behaviour of c(F) looks quite smooth, this giving evidence for 
the physics phenomenon we are studying being not a crystallization process. Moreover, we 
never observed the system definitely trapped in a configuration as it happens in the previously 
discussed N=32 case. 

We plot in [Fig. 11] data on P{d) for different A^ values at the highest temperature 
considered (F=l). The shape is the Gaussian-like one characteristic of the liquid phase and 
it changes very little when varying the number of particles. 

Finally we present in [Figs. 12-15] the low temperature data on P{d) for A^=28, 30, 34 
and 36, respectively. The qualitative features are the same for different number of particles. 
The behaviour changes in a remarkably short F range around F = F*, and P{d) looks highly 
non trivial at higher F values. The physical meaning of the peaks shown by the equilibrium 
probability distribution of distance between states in the low temperature region comes from 
the mean field theoretical picture. In this region the data are well consistent with the scenario 
that in the glassy phase a small number of valleys in the free energy landscape give the 
dominant contribution to the partition function. 

Both the value of F* and the P{d) shape in the glassy phase are strongly dependent on 
N. Differences are evident also when the number of particles is varied only of 2 (i.e. between 
N=28 and 30 or A^=34 and 36). This is easily understandable since a little difference in A^ 
can change abruptly the kinds (and the number) of configurations that maximize and that are 
near to maximize relative distances between particles (minimizing the Hamiltonian) . On the 
other hand, the observed non trivial dependence of the behaviour of the model on the number 
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Figure 10: Data on the specific heath c as a function of F for A^=28(+), 30(x), 34(*) and 
36(n). Lines are only to join neighbouring points. 
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Figure 11: Data on P{d) at F=l for A^=28(+), 30(x), 34(>k) and 36(n) 
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Figure 12: Data on P{d) for A^=28. From left to right and top to bottom r=1.5, 1.6, 1.7, 1.8, 
1.9 and 2. Plotted data have been measured during a 2'^^PT steps run but we stress that they 
are perfectly compatible with those obtained by a previous 2^^ step run. 
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Figure 13: Data on P{d) for A^=30. From left to right and top to bottom r=1.5, 1.6, 1.7, 1.8, 
1.9 and 2 (i. e. the same value that in previous [Fig. 10]). Here data have been measured 
during a 2^^PT step run. 
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Figure 14: Data on P{d) for N=34. From left to right and top to bottom r=1.2, 1.3, 1.4, 1.6, 
1.8 and 2. Here and in the next case of N=36 data have been measured during a 2^'^PT step 



run. 
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Figure 15: Data on P{d) for A^=36. From left to right and top to bottom r=1.2, 1.3, 1.4, 1.6, 
1.8 and 2.0 (i. e. the same value that in previous [Fig. 15]). 
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of particles seems to agree with A^ playing in some way in structural glasses the same role that 
quenched random variables do in spin glasses. 

As a last remark we note that when reaching F values well higher than F* some peaks of 
P{d) disappear and the height of the other ones (usually corresponding to smaller d values) 
goes up. This seems quite reasonable in a finite size system with continuous degrees of freedom 
since in spite of the absence of a perfectly crystalline ground state only a few configurations 
of the particles are expected to really minimize the Hamiltonian. 

The natural next step in this analysis should consist in trying to get the behaviour of 



quantities averaged over different A^ values. We advance that from early new results p9[ a 
slightly different definition of distance seems to be more suitable for mixtures and it could make 
easier this kind of study, permitting to gain a further insight into structural glass properties. 
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